Hands-on Exercise 3

Author

Xu Haiyang

Published

August 26, 2024

Modified

August 28, 2024

1. Overview

Spatial Point Pattern Analysis examines the distribution pattern of points on a surface. These points can represent the locations of:

  • Events, such as crimes, traffic accidents, or disease outbreaks.

  • Business services, like coffee shops and fast food outlets, or facilities such as childcare and eldercare centers.

In this exercise, I utilize functions from the spatstat package to explore the spatial point processes of childcare centers in Singapore.

Key Questions:

  • Are the childcare centers in Singapore randomly distributed across the country?

  • If not, where are the areas with a higher concentration of childcare centers?

2. The Data

To answer these questions, I used three datasets:

  • CHILDCARE: A point dataset with location and attribute information for childcare centers in Singapore. This dataset, in GeoJSON format, was downloaded from Data.gov.sg.

  • MP14_SUBZONE_WEB_PL: A polygon dataset representing the 2014 Master Plan Planning Subzone boundaries provided by the Urban Redevelopment Authority (URA) in ESRI Shapefile format, also from Data.gov.sg.

  • CostalOutline: A polygon dataset outlining the national boundary of Singapore, provided by the Singapore Land Authority (SLA) in ESRI Shapefile format.

3. Installing and Loading the R Packages

I installed and loaded the necessary R packages to handle spatial data, perform point pattern analysis, and create thematic maps:

pacman::p_load(sf, raster, spatstat, tmap, tidyverse)
  • pacman::p_load(): Ensures that the required packages are installed and loaded. The packages include:

    • sf: For handling and analyzing spatial vector data.

    • raster: For raster data manipulation.

    • spatstat: For spatial point pattern analysis.

    • tmap: For creating thematic maps.

    • tidyverse: For general data manipulation and visualization.

4. Spatial Data Wrangling

4.1 Importing the Spatial Data

I imported the datasets and ensured they all use the same coordinate reference system (CRS) for consistency in spatial analysis.

childcare_sf <- st_read("data/child-care-services-geojson.geojson") %>%
  st_transform(crs = 3414)
Reading layer `child-care-services-geojson' from data source 
  `C:\EasonXu-HY99\IS415\Hands-on_Ex\Hands-on_Ex03\data\child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
  • st_transform(): Converts the spatial data to the Singapore-specific projected CRS (EPSG: 3414).
sg_sf <- st_read(dsn = "data", layer = "CostalOutline") 
Reading layer `CostalOutline' from data source 
  `C:\EasonXu-HY99\IS415\Hands-on_Ex\Hands-on_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21
mpsz_sf <- st_read(dsn = "data", layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\EasonXu-HY99\IS415\Hands-on_Ex\Hands-on_Ex03\data' using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
st_crs(childcare_sf)  
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(mpsz_sf)       
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(sg_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
  • st_crs(): Checks the CRS of each dataset to ensure they are compatible for spatial operations.
mpsz_sf <- st_set_crs(mpsz_sf, 3414)
Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that
sg_sf <- st_set_crs(sg_sf, 3414)
Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
that
  • st_set_crs(): Assigns the specified CRS (EPSG: 3414) to the datasets if not already defined, ensuring all datasets are in the same spatial reference.
childcare_sf <- st_transform(childcare_sf, crs = 3414)
  • Re-applied the transformation to make sure the childcare dataset uses the correct CRS.
st_crs(childcare_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
st_crs(sg_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
  • Confirmed the CRS settings for all datasets.

4.2 Mapping the Geospatial Data Sets

To visualize the spatial datasets, I created static and interactive maps using the tmap package:

tmap_mode("plot")
tmap mode set to plotting
tm_shape(sg_sf) +
  tm_polygons(col = "grey", border.col = "black") +
tm_shape(mpsz_sf) +  
  tm_polygons(col = "grey", border.col = "black") +
tm_shape(childcare_sf) +  
  tm_dots(col = "black", size = 0.1)

  • tmap_mode("plot"): Sets the mode for static plotting.

  • tm_shape(): Specifies the spatial object to be used in the map.

  • tm_polygons(): Plots polygon features with specified fill and border colors.

  • tm_dots(): Adds point symbols to the map for representing the childcare centers.

For an interactive map:

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare_sf) + tm_dots()
  • tmap_mode('view'): Switches to interactive viewing mode, allowing dynamic exploration of the spatial data.

Returning to static plotting mode:

tmap_mode('plot')
tmap mode set to plotting
  • Resets to static plotting after using the interactive mode.

5. Geospatial Data Wrangling

5.1 Converting sf Data Frames to sp’s Spatial* Class

To perform certain spatial analyses, I needed to convert sf data frames into sp’s Spatial* class objects.

childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)
  • as_Spatial(): Converts an sf object into an sp object.
summary(childcare)
Object of class SpatialPointsDataFrame
Coordinates:
               min      max
coords.x1 11203.01 45404.24
coords.x2 25667.60 49300.88
coords.x3     0.00     0.00
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
+x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs]
Number of points: 1545
Data attributes:
     Name           Description       
 Length:1545        Length:1545       
 Class :character   Class :character  
 Mode  :character   Mode  :character  
summary(mpsz)
Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x  2667.538 56396.44
y 15748.721 50256.33
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
+x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs]
Data attributes:
    OBJECTID       SUBZONE_NO      SUBZONE_N          SUBZONE_C        
 Min.   :  1.0   Min.   : 1.000   Length:323         Length:323        
 1st Qu.: 81.5   1st Qu.: 2.000   Class :character   Class :character  
 Median :162.0   Median : 4.000   Mode  :character   Mode  :character  
 Mean   :162.0   Mean   : 4.625                                        
 3rd Qu.:242.5   3rd Qu.: 6.500                                        
 Max.   :323.0   Max.   :17.000                                        
    CA_IND           PLN_AREA_N         PLN_AREA_C          REGION_N        
 Length:323         Length:323         Length:323         Length:323        
 Class :character   Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character   Mode  :character  
                                                                            
                                                                            
                                                                            
   REGION_C           INC_CRC            FMEL_UPD_D             X_ADDR     
 Length:323         Length:323         Min.   :2014-12-05   Min.   : 5093  
 Class :character   Class :character   1st Qu.:2014-12-05   1st Qu.:21864  
 Mode  :character   Mode  :character   Median :2014-12-05   Median :28465  
                                       Mean   :2014-12-05   Mean   :27257  
                                       3rd Qu.:2014-12-05   3rd Qu.:31674  
                                       Max.   :2014-12-05   Max.   :50425  
     Y_ADDR        SHAPE_Leng        SHAPE_Area      
 Min.   :19579   Min.   :  871.5   Min.   :   39438  
 1st Qu.:31776   1st Qu.: 3709.6   1st Qu.:  628261  
 Median :35113   Median : 5211.9   Median : 1229894  
 Mean   :36106   Mean   : 6524.4   Mean   : 2420882  
 3rd Qu.:39869   3rd Qu.: 6942.6   3rd Qu.: 2106483  
 Max.   :49553   Max.   :68083.9   Max.   :69748299  
summary(sg)
Object of class SpatialPolygonsDataFrame
Coordinates:
        min      max
x  2663.926 56047.79
y 16357.981 50244.03
Is projected: TRUE 
proj4string :
[+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
+x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs]
Data attributes:
    GDO_GID          MSLINK          MAPID    COSTAL_NAM       
 Min.   : 1.00   Min.   : 1.00   Min.   :0   Length:60         
 1st Qu.:15.75   1st Qu.:17.75   1st Qu.:0   Class :character  
 Median :30.50   Median :33.50   Median :0   Mode  :character  
 Mean   :30.50   Mean   :33.77   Mean   :0                     
 3rd Qu.:45.25   3rd Qu.:49.25   3rd Qu.:0                     
 Max.   :60.00   Max.   :67.00   Max.   :0                     

5.2 Converting the Spatial* Class into Generic sp Format

I further converted the sp objects into more generic formats used by other spatial analysis functions.

childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")
  • Why: Converting to generic sp formats allows compatibility with a wider range of functions and analyses in R.

  • Functions:

    • as(): Converts objects from one class to another, in this case, converting to “SpatialPoints” and “SpatialPolygons”.

Checking the converted objects:

childcare_sp
class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sg_sp
class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
  • Challenge: Understanding the differences between the Spatial* classes and generic sp objects is crucial for using the correct format in different analyses.

5.3 Converting the Generic sp Format into spatstat’s ppp Format

To analyze spatial point patterns using spatstat, I converted the generic sp objects into ppp (planar point pattern) format.

childcare_ppp <- as.ppp(childcare_sf)
Warning in as.ppp.sf(childcare_sf): only first attribute column is used for
marks
childcare_ppp
Marked planar point pattern: 1545 points
marks are of storage type  'character'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
  • Why: The ppp format is specifically designed for spatial point pattern analysis in spatstat.

  • Functions:

    • as.ppp(): Converts spatial objects into the ppp format.

Visualizing and summarizing the ppp object:

plot(childcare_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

summary(childcare_ppp)
Marked planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units
  • plot(): Visualizes the spatial distribution of points in the ppp object.

  • summary(): Provides a detailed summary of the ppp object, including the number of points and window properties.

5.4 Handling Duplicated Points

I checked for and handled any duplicated points in the dataset to ensure the accuracy of the spatial analysis.

any(duplicated(childcare_ppp))
[1] FALSE
  • Why: Detecting duplicated points is crucial as they can skew spatial point pattern analyses.

  • Functions:

    • any(): Checks if there are any TRUE values in a logical vector, indicating duplicates in this case.

    • duplicated(): Identifies duplicated points in the ppp object.

Checking the multiplicity of points:

multiplicity(childcare_ppp)
   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1037] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1074] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1111] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1148] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1185] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1222] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1259] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1296] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1333] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1370] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1407] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1444] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1481] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1518] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
sum(multiplicity(childcare_ppp) > 1)
[1] 0
  • multiplicity(): Returns the number of times each point occurs.

  • sum(): Sums the total number of duplicate points.

Visualizing duplicated points:

tmap_mode('view')
tmap mode set to interactive viewing
tm_shape(childcare) + tm_dots(alpha=0.4, size=0.05)
tmap_mode('plot')
tmap mode set to plotting
  • Why: Visualizing the data helps identify and understand the location and extent of duplicated points.

Handling duplicates by jittering points:

childcare_ppp_jit <- rjitter(childcare_ppp, retry=TRUE, nsim=1, drop=TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE
  • rjitter(): Randomly displaces points to reduce overlap, helping to handle duplicates while retaining the general spatial pattern.

5.5 Creating owin Object

I created an owin object to define the observation window for point pattern analysis.

sg_owin <- as.owin(sg_sf)
  • Why: Defining an observation window is necessary for controlling the area within which spatial point patterns are analyzed.

  • Functions:

    • as.owin(): Converts a spatial object into an owin object, defining a spatial observation window.

Plotting and summarizing the owin object:

plot(sg_owin)

summary(sg_owin)
Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401

5.6 Combining Point Events Object and owin Object

I combined the ppp and owin objects to focus on the analysis within the specified spatial boundary.

childcareSG_ppp = childcare_ppp[sg_owin]
summary(childcareSG_ppp)
Marked planar point pattern:  1545 points
Average intensity 2.129929e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401
plot(childcareSG_ppp)
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 1545 symbols are shown in the symbol map

  • Why: Combining the point events with the observation window ensures that the analysis is conducted only within the desired spatial extent.

6. First-Order Spatial Point Patterns Analysis

6.1 Kernel Density Estimation

6.1.1 Computing Kernel Density Estimation Using Automatic Bandwidth Selection Method
kde_childcareSG_bw <- density(childcareSG_ppp, sigma=bw.diggle, edge=TRUE, kernel="gaussian") 
plot(kde_childcareSG_bw)

bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
298.4095 
  • Why: Kernel Density Estimation (KDE) helps identify areas with higher point concentration. The bandwidth selection method (bw.diggle) is used to optimize the KDE for the data’s distribution.

  • Functions:

    • density(): Computes the KDE for the point pattern.

    • bw.diggle(): Selects the bandwidth automatically using Diggle’s method, which is suitable for spatial data with edge corrections.

6.1.2 Rescaling KDE Values
childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")
kde_childcareSG.bw <- density(childcareSG_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG.bw)

  • Why: Rescaling helps in interpreting the KDE in a more familiar unit (e.g., kilometers).

  • Functions:

    • rescale.ppp(): Rescales the spatial coordinates of the ppp object.

6.2 Working with Different Automatic Bandwidth Methods

bw.CvL(childcareSG_ppp.km)
   sigma 
4.543278 
bw.scott(childcareSG_ppp.km)
 sigma.x  sigma.y 
2.224898 1.450966 
bw.ppl(childcareSG_ppp.km)
    sigma 
0.3897114 
bw.diggle(childcareSG_ppp.km)
    sigma 
0.2984095 
  • Why: Exploring different bandwidth methods allows for comparison and selection of the best fit for the data.

  • Functions:

    • bw.CvL(), bw.scott(), bw.ppl(), bw.diggle(): Different methods for selecting the bandwidth in KDE.

Visualizing KDE with different bandwidths:

kde_childcareSG.ppl <- density(childcareSG_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

6.3 Working with Different Kernel Methods

par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km, sigma=bw.ppl, edge=TRUE,  kernel="gaussian"), main="Gaussian")
plot(density(childcareSG_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="epanechnikov"), main="Epanechnikov")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="quartic"), main="Quartic")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
plot(density(childcareSG_ppp.km, sigma=bw.ppl, edge=TRUE, kernel="disc"), main="Disc")
Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel

  • Why: Different kernel functions can affect the KDE results. Comparing different kernels helps determine the most appropriate method for the data.

  • Functions:

    • density() with different kernel options (gaussian, epanechnikov, quartic, disc): Computes KDE using different smoothing kernels.

7. Fixed and Adaptive KDE

7.1 Computing KDE Using Fixed Bandwidth

kde_childcareSG_600 <- density(childcareSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600)

  • Why: A fixed bandwidth KDE uses a constant smoothing parameter (sigma) across the entire study area, providing a uniform level of smoothing. This approach is useful for identifying general patterns in the distribution of points.

  • Functions:

    • density(): Computes kernel density estimates for point patterns. The sigma parameter specifies the bandwidth.

7.2 Computing KDE Using Adaptive Bandwidth

kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)

  • Why: An adaptive bandwidth KDE adjusts the smoothing parameter based on local point density, providing finer detail in areas with high point concentration and smoother estimates in sparse areas. This method is useful for identifying local clusters.

  • Functions:

    • adaptive.density(): Computes KDE with adaptive bandwidth, adjusting the bandwidth according to point density.

Comparing fixed and adaptive bandwidth:

par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")

  • Why: Comparing both methods visually helps to understand how different bandwidth strategies affect the KDE.

7.3 Converting KDE Output into Grid Object

gridded_kde_childcareSG_bw <- as(kde_childcareSG.bw, "SpatialGridDataFrame")
spplot(gridded_kde_childcareSG_bw)

  • Why: Converting KDE output into a grid format allows for easier visualization and manipulation within various GIS tools.

  • Functions:

    • as(): Converts objects from one class to another.

    • spplot(): Creates spatial plots of gridded data.

7.3.1 Converting Gridded Output into Raster
kde_childcareSG_bw_raster <- raster(kde_childcareSG.bw)
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -8.476185e-15, 28.51831  (min, max)
  • Why: Raster conversion facilitates the integration of KDE results with other raster-based analyses or visualization techniques.

  • Functions:

    • raster(): Converts a spatial object into a raster format.
7.3.2 Assigning Projection Systems
projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : layer 
values     : -8.476185e-15, 28.51831  (min, max)
  • Why: Assigning the correct CRS ensures that spatial analyses and visualizations are accurate and geographically meaningful.

  • Functions:

    • projection(): Assigns or retrieves the CRS of a raster object.

7.4 Visualizing the Output in tmap

tm_shape(kde_childcareSG_bw_raster) + 
  tm_raster("layer", palette = "viridis") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

  • Why: Visualizing KDE results with tmap provides a more intuitive and aesthetically pleasing representation of spatial data.

  • Functions:

    • tm_shape(): Specifies the spatial object to be visualized.

    • tm_raster(): Visualizes raster data with a color gradient.

    • tm_layout(): Customizes the layout of the map.

7.5 Comparing Spatial Point Patterns Using KDE

7.5.1 Extracting Study Areas
pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")
  • Why: Extracting specific study areas allows focused analysis on different regions to compare spatial point patterns.

  • Functions:

    • filter(): Subsets data based on specified conditions.

Plotting study areas:

par(mfrow=c(2,2))
plot(pg, main = "Ponggol")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(tm, main = "Tampines")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

plot(ck, main = "Choa Chu Kang")
Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all

plot(jw, main = "Jurong West")
Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all

7.5.2 Creating owin Object
pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)
  • Why: Converting study areas into owin objects defines the observation window for point pattern analysis within each area.

  • Functions:

    • as.owin(): Converts a spatial object into an owin format.
7.5.3 Combining Childcare Points and the Study Area
childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]
  • Why: Combining points with their respective study areas allows for localized point pattern analysis.

Rescaling point patterns:

childcare_pg_ppp.km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale.ppp(childcare_jw_ppp, 1000, "km")
  • Why: Rescaling the data to kilometers facilitates comparison across different study areas.

  • Functions:

    • rescale.ppp(): Rescales the coordinates of point pattern objects.

Visualizing rescaled point patterns:

par(mfrow=c(2,2))
plot(childcare_pg_ppp.km, main="Punggol")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_tm_ppp.km, main="Tampines")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 89 symbols are shown in the symbol map
plot(childcare_ck_ppp.km, main="Choa Chu Kang")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 61 symbols are shown in the symbol map
plot(childcare_jw_ppp.km, main="Jurong West")
Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different character
Warning: Only 10 out of 88 symbols are shown in the symbol map

7.5.4 Computing KDE
par(mfrow=c(2,2))
plot(density(childcare_pg_ppp.km, sigma=bw.diggle, edge=TRUE,  kernel="gaussian"), main="Punggol")
plot(density(childcare_tm_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian"), main="Tempines")
plot(density(childcare_ck_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian"), main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian"), main="Jurong West")

  • Why: Computing KDE for each area helps identify and compare spatial distribution patterns within different regions.
7.5.5 Computing Fixed Bandwidth KDE
par(mfrow=c(2,2))
plot(density(childcare_ck_ppp.km, sigma=0.25, edge=TRUE,  kernel="gaussian"), main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km, sigma=0.25, edge=TRUE, kernel="gaussian"), main="Jurong West")
plot(density(childcare_pg_ppp.km, sigma=0.25, edge=TRUE, kernel="gaussian"), main="Punggol")
plot(density(childcare_tm_ppp.km, sigma=0.25, edge=TRUE,  kernel="gaussian"), main="Tampines")

  • Why: Using a fixed bandwidth allows for direct comparison of KDE results across different study areas with a uniform smoothing parameter.

8. Nearest Neighbor Analysis

8.1 Testing Spatial Point Patterns Using Clark and Evans Test

clarkevans.test(childcareSG_ppp, correction="none", clipregion="sg_owin", alternative=c("clustered"), nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.55631, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
  • Why: The Clark and Evans test determines whether a point pattern is more clustered, random, or regular compared to a Poisson distribution.

  • Functions:

    • clarkevans.test(): Performs the Clark and Evans test for spatial point patterns.

8.2 Clark and Evans Test: Choa Chu Kang Area

clarkevans.test(childcare_ck_ppp, correction="none", clipregion=NULL, alternative=c("two.sided"), nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.95577, p-value = 0.5087
alternative hypothesis: two-sided

8.3 Clark and Evans Test: Tampines Planning Area

clarkevans.test(childcare_tm_ppp, correction="none", clipregion=NULL, alternative=c("two.sided"), nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.77884, p-value = 6.569e-05
alternative hypothesis: two-sided
  • Why: Performing the Clark and Evans test on different areas allows for localized analysis of spatial point patterns to identify variations in clustering.

Continue to Page 2